7393 Exam 2 R Code Andrei Matveev

library(data.table)
library(truncnorm)
library(coda)
devtools::install_github("stan-dev/loo")
library(loo)
library(rstan) 
library(shinystan)
library(invgamma)
library(bayesplot)
library(posterior)

devtools::install_github("jfrench/bayesutils")

setwd('/Users/AM/Documents/_CU Masters/2020 fall Bayesian_7393/code/Bayesian_Statistics_Class_Code/Exam_code')
#df1 = load("diamonds_simple.rda")
#rm(df1, diamonds_simple)
df = bayesutils::diamonds_simple

df$lprice = log(df$price)
df$lcarat = log(df$carat)

lq_theta_y = function(sigmaSQ, beta0, beta1, lpr = df$lprice, lcar = df$lcarat) {
  ld = dnorm(beta0, 0, 100, log = TRUE) + 
    dnorm(beta1, 0, 100, log = TRUE) + 
    dinvgamma(sigmaSQ, shape = 0.01, rate = 0.01, log = TRUE) +
    sum(dnorm(lpr, mean = (beta0 + beta1 * lcar), sd = sqrt(sigmaSQ), log = TRUE))
  return(ld)
}

mh = function(B, theta_start) {
  theta = array(0, c((B+1), 3), dimnames = list(c(), c("sigmaSQ", "beta0", "beta1")))
  
  theta[1,1] = theta_start[1]
  theta[1,2] = theta_start[2]
  theta[1,3] = theta_start[3]

  for (i in 2:dim(theta)[1]) {
    
    ###  step for sigmaSQ
    beta0_star = theta[(i-1),2]
    beta1_star = theta[(i-1),3]
    
    sigmaSQ_star = rtruncnorm(n = 1, a=0, b=Inf, mean = theta[(i-1),1], sd = 0.1)

    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      log(dtruncnorm(x = sigmaSQ_star, a=0, b=Inf, mean = theta[(i-1),1], sd = 0.1))
    den_logr = lq_theta_y(sigmaSQ = theta[i-1, 1], beta0 = beta0_star, beta1 = beta1_star) -
      log(dtruncnorm(x = theta[i-1, 1], a=0, b=Inf, mean = sigmaSQ_star, sd = 0.1)) 
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,1] = sigmaSQ_star
    } else {
      theta[i,1] = theta[(i - 1), 1]
    }
    
    ###  step for beta0
    beta1_star = theta[(i-1),3] # it is the repeated code, but I need it to keep the interpretability
    sigmaSQ_star = theta[i,1] # update sigmaSQ_star after the Gibbs step for sigmaSQ
    
    beta0_star = rnorm(1, theta[(i-1),2], 0.1) # !!! check the parametrization (0.1 or 0.1^2)
    
    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      dnorm(x = beta0_star, mean = theta[(i-1),2], sd = 0.1, log = TRUE)
    den_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = theta[i-1, 2], beta1 = beta1_star) -
      dnorm(x = theta[(i-1),2], mean = beta0_star, sd = 0.1, log = TRUE)
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,2] = beta0_star
    } else {
      theta[i,2] = theta[(i - 1), 2]
    }
    
    ###  step for beta1
    sigmaSQ_star = theta[i,1] # it is the repeated code, but I need it to keep the interpretability
    beta0_star = theta[i,2] # update beta0_star after the Gibbs step for beta0

    beta1_star = rnorm(1, theta[(i-1),3], 0.1) # !!! check the parametrization (0.1 or 0.1^2)
    
    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      dnorm(x = beta1_star, mean = theta[(i-1),3], sd = 0.1, log = TRUE)
    den_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = theta[i-1, 3]) -
      dnorm(x = theta[(i-1),3], mean = beta1_star, sd = 0.1, log = TRUE)
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,3] = beta1_star
    } else {
      theta[i,3] = theta[(i - 1), 3]
    }
    
  }
  return(theta)
}

B = 10^5 
keep = (B/2 + 1):(B + 1)
chain1 = mh(B, theta_start = c(0.1, -1, -1))
chain2 = mh(B, theta_start = c(0.3, 0, 0))
chain3 = mh(B, theta_start = c(0.5, -1, 1))
chain4 = mh(B, theta_start = c(0.2, 1, 1))

mc = mcmc.list(mcmc(chain1[keep,]), mcmc(chain2[keep,]),
                 mcmc(chain3[keep,]), mcmc(chain4[keep,]))
summary(mc)

Iterations = 1:50001
Thinning interval = 1 
Number of chains = 4 
Sample size per chain = 50001 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean      SD  Naive SE Time-series SE
sigmaSQ 0.1152 0.00925 2.068e-05      6.977e-05
beta0   8.3798 0.02240 5.008e-05      1.647e-04
beta1   1.5064 0.03135 7.011e-05      2.141e-04

2. Quantiles for each variable:

           2.5%    25%    50%    75%  97.5%
sigmaSQ 0.09853 0.1088 0.1145 0.1211 0.1348
beta0   8.33546 8.3649 8.3797 8.3948 8.4234
beta1   1.44522 1.4854 1.5064 1.5276 1.5680

gelman.diag(mc, autoburnin = FALSE)

Potential scale reduction factors:

        Point est. Upper C.I.
sigmaSQ       1.03       1.08
beta0         1.00       1.01
beta1         1.01       1.03

Multivariate psrf

1.03

# extracting the MCMC samples from the soda_fit object
samples = extract(r_fit_A)
ncycles = length(samples[[1]])
T <- length(df$price)

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$lcarat[i]
  yrep[, i] = rnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

# approximate posterior predictive density for an observation
# with the same covariates as observation 1
#plot(density(yrep[,1]))
color_scheme_set("red")

# posterior predictive check
y = df$lprice
ppc_hist(y, yrep[sample(1:ncycles, 8),]) + ggtitle("Model A: histogram comparing the response to 8 replicated data sets")

# scatterplot of y vs average yrep
ppc_scatter_avg(y, yrep) + ggtitle("Model A: scatterplot comparing the response to the average replicated data set")

ppc_dens_overlay(y, yrep = yrep[sample(1:ncycles, 100),]) + ggtitle("Model A: density plot comparing the density of the response to the densities of 100 replicated data sets")

ppc_error_scatter_avg(y, yrep = yrep) + ggtitle("Model A: scatter plot comparing the response to the average predictive error")

ppc_intervals(y, yrep = yrep)

ppc_error_scatter_avg_vs_x(y = y, yrep = yrep, x = df$lcarat)

ppc_intervals(y = y, yrep = yrep, x = df$lcarat) + ggtitle("Model A: Posterior predictive intervals to the observed data values") + xlab("log carat") + ylab("log price")

# extracting the MCMC samples from the soda_fit object
samples = extract(r_fit_B)
ncycles = length(samples[[1]])
T <- length(df$price)

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$carat[i] + samples$beta2 * df$carat[i] * df$carat[i]
  yrep[, i] = rnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

# approximate posterior predictive density for an observation
# with the same covariates as observation 1
#plot(density(yrep[,1]))
color_scheme_set("green")

# posterior predictive check
y = df$price
ppc_hist(y, yrep[sample(1:ncycles, 8),]) + ggtitle("Model B: histogram comparing the response to 8 replicated data sets")

# scatterplot of y vs average yrep
ppc_scatter_avg(y, yrep) + ggtitle("Model B: scatterplot comparing the response to the average replicated data set")

ppc_dens_overlay(y, yrep = yrep[sample(1:ncycles, 100),]) + ggtitle("Model B: density plot comparing the density of the response to the densities of 100 replicated data sets")

ppc_error_scatter_avg(y, yrep = yrep) + ggtitle("Model B: scatter plot comparing the response to the average predictive error")

ppc_intervals(y = y, yrep = yrep, x = df$carat) + ggtitle("Model B: Posterior predictive intervals to the observed data values") + xlab("carat") + ylab("price")